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Abstract 

This  paper  describes  a  simple  isothermal  two-phase  flow  dynamic  model  that  predicts  the  experimentally  observed  temporal  behavior  of  a  proton 
exchange  membrane  fuel  cell  stack.  This  model  is  intended  for  use  in  embedded  real  time  control  where  computational  simplicity  is  of  critical 
importance.  A  reproducible  methodology  is  presented  to  experimentally  identify  six  (6)  tunable  physical  parameters  based  on  the  estimation  of 
the  cell  voltage,  the  water  vapor  transport  through  the  membrane  and  the  accumulation  of  liquid  water  in  the  gas  channels.  The  model  equations 
allow  temporal  calculation  of  the  species  concentrations  across  the  gas  diffusion  layers,  the  vapor  transport  across  the  membrane,  and  the  degree 
of  flooding  within  the  cell  structure.  The  notion  of  apparent  current  density  then  relates  this  flooding  phenomena  to  cell  performance  through  a 
reduction  in  the  cell  active  area  as  liquid  water  accumulates.  Despite  the  oversimplification  of  many  complex  phenomena,  this  model  provides  a 
useful  tool  for  predicting  the  resulting  decay  in  cell  voltage  over  time  only  after  it  has  been  tuned  with  experimental  data.  The  calibrated  model 
and  tuning  procedure  is  demonstrated  with  a  1.4  kW  (24  cell,  300  cm2)  stack,  using  pressure  regulated  pure  hydrogen  supplied  to  a  dead-ended 
anode,  under  a  range  of  operating  conditions  typical  for  multi-cell  stacks. 

©  2007  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

The  management  of  water  is  critical  for  optimizing  perfor¬ 
mance  of  a  polymer  electrolyte  membrane  fuel  cell  (PEMFC) 
stack.  Because  the  ionic  conductivity  of  the  membrane  is  depen¬ 
dent  upon  its  water  content  [1],  a  balance  must  be  struck 
between  reactant  (hydrogen  and  oxygen)  delivery  and  water 
supply  and  removal.  Depending  upon  the  operating  conditions 
of  the  PEMFC  stack,  the  flow  patterns  in  the  anode  and  cath¬ 
ode  channels,  and  the  design  of  the  anode  gas  delivery  system 
(dead-ended  or  flow  through),  this  liquid  water  can  accumulate 
within  the  gas  diffusion  layers  (GDLs)  and  channels  [2-5],  as 
shown  in  Fig.  1.  Whether  obstructing  reactant  flow  or  reduc¬ 
ing  the  number  of  active  catalyst  sites,  the  impact  of  flooding 


*  Corresponding  author.  Tel.:  +1  7347644272. 

E-mail  address:  dmckay@umich.edu  (D.A.  McKay). 

1  Currently  with  the  Southwest  Research  Institute,  USA. 

0378-7753/$  -  see  front  matter  ©  2007  Elsevier  B.V.  All  rights  reserved, 
doi:  10. 1016/j.jpowsour.2007. 12.031 


is  a  reduction  in  the  power  output  of  the  fuel  cell  stack,  seen 
by  a  decrease  in  cell  voltage  [6].  Thus,  a  real-time  estima¬ 
tion  of  the  degree  of  flooding  within  the  cell  structure  and  its 
impact  on  the  cell  electrical  output  with  standard,  cheap,  and 
reliable  sensors  is  critical  for  active  water  management.  More¬ 
over,  a  low  order  control-oriented  model  must  be  derived  for 
further  considering  such  issues  as  identifiability,  observability, 
and  controllability. 

To  gain  a  better  understanding  of  reactant  and  water  transport 
within  the  GDL  and  catalyst  layers,  many  CFD  models  have 
been  developed  to  approximate  the  two-  or  three-dimensional 
flow  of  hydrogen,  air,  and  water  at  steady-state  within  the 
cell  structure  [7-12].  Using  experimental  steady-state  polariza¬ 
tion  (voltage  versus  current)  data  for  parameter  identification, 
Guo  et  al.  [13]  and  Carnes  and  Djilali  [14]  investigated  the 
sensitivity  of  the  cell  performance  to  the  identified  param¬ 
eters.  Further,  using  a  model  to  simulate  polarization  data 
with  a  given  set  of  parameters,  constrained  quadratic  program¬ 
ming  was  then  used  to  identify  these  given  parameters  [15] 
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Fig.  1,  Schematic  of  capillary  flow  of  liquid  water  through  the  gas  diffusion 
layers. 

and  address  parameter  identifiability  and  uniqueness  issues 
[16]. 

While  these  models  are  ideal  for  investigating  transport  phe¬ 
nomena  with  two  phase  flow  and  spatial  gradients,  examining 
parameter  sensitivity,  or  the  influence  of  material  properties 
on  cell  performance,  experimental  validation  of  these  models, 
often  completed  by  comparing  measured  to  estimated  polariza¬ 
tion  curves,  is  still  lacking.  A  few  publications  with  steady- state 
validation  efforts  (i)  point  to  a  mismatch  between  model  predic¬ 
tion  and  spatially  resolved  experimental  data  [17]  indicating  that 
different  spatial  distributions  can  correspond  to  a  similar  aver¬ 
aged  polarization  curve  [16,17],  and  (ii)  achieve  good  prediction 
of  steady-state  and  spatially  resolved  current  density  measure¬ 
ments  after  tuning  parameters  to  several  orders  of  magnitude  of 
their  theoretical  values  [18]. 

Although  steady-state  polarization  measurements  do  not 
offer  a  conclusive  data  set  for  model  validation,  the  transient 
polarization  response  provides  useful  data  for  model  valida¬ 
tion  especially  during  unsteady  operation  such  as  flooding 
[19,20].  Several  transient  models  have  been  reported  to  illicit 
the  relationship  between  critical  material  properties  and  oper¬ 
ating  conditions  on  the  dynamic  fuel  cell  response  [21-25], 
however  few  have  been  validated  against  transient  experimental 
data  and  are  of  sufficient  complexity  for  implementation  in  real 
time  control  applications. 

Control-oriented  transient  models  have  been  developed  to 
account  for  the  formation  of  liquid  water  within  the  gas  chan¬ 
nels  [26]  or  within  both  the  channels  and  the  GDL  [27],  however 
they  do  not  relate  the  effect  of  flooding  to  decreased  cell  poten¬ 
tial,  a  key  indication  of  how  flooding  impacts  cell  performance. 
A  relationship  between  flooding  and  cell  performance  was  intro¬ 
duced  in  Ref.  [20],  appearing  later  in  Ref.  [28],  using  the  notion 
of  apparent  current  density  to  relate  the  accumulation  of  liquid 
water  in  the  gas  channels  to  a  reduction  in  the  cell  active  area, 
in  turn  increasing  the  cell  current  density  and  lowering  cell  volt¬ 
age.  Although  the  apparent  current  density  calculation  based  on 
the  water  accumulation  in  the  channels  approximates  the  cell 
voltage  behavior  well  during  a  range  of  transient  and  steady 
conditions  the  stack  typically  operates  in,  more  experimental 
evidence  and  justification  of  this  simplification  is  needed  and  is 
underway  in  our  laboratory. 


In  this  paper,  we  present  a  low  order  model  of  the  liquid 
water  and  gas  dynamics  within  the  GDL  to  simulate  both  the 
effects  of  reactant  starvation  and  flooding.  We  focus  on  the  one¬ 
dimensional  dynamics  through  the  GDL  thickness,  assuming 
invariant  two-dimensional  properties  in  the  plane  parallel  to 
the  membrane,  lumped  volume  manifold  filling  dynamics  for 
the  gas  channels,  and  lumped  parameter  characteristics  for  the 
membrane.  Lumping  the  GDL  and  channel  into  a  single  vol¬ 
ume,  Hernandez  et  al.  [29]  experimentally  validated  their  model 
for  a  flow-through  anode  with  no  gas  dynamics.  Lumping  the 
GDL  volume  was  also  pursued  in  Ref.  [28]  and  validated  against 
experimental  data  for  a  Ballard®  NEXA  ™  system.  Note  that 
[28]  nearly  doubled  the  number  of  experimentally  identified 
parameters  from  the  work  originally  presented  in  Ref.  [20]  and 
used  here. 

In  this  paper  we  extend  and  test  the  validity  of  [20]  to  a 
wider  range  of  current  densities  (0-0.3  A  cm-2),  temperature 
(45-65  °C)  and  air  stoichiometries  (150-300%).  These  con¬ 
ditions  are  tested  while  the  stack  operates  mostly  under  full 
hydrogen  utilization  with  intermittent  and  short  high  hydrogen 
flow  conditions  associated  with  dead-ended  anode  operation.  It 
is  shown  that  our  model  predicts  both  the  fast  voltage  dynamics 
during  step  changes  in  current  (the  gas  dynamics)  and  the  slow 
voltage  behavior  while  liquid  water  is  accumulating  in  the  GDL 
and  gas  channels  (water  dynamics),  whereas  [28]  only  predicts 
the  slow  voltage  dynamics  well.  Hence,  the  model  presented 
here  can  be  used  for  estimation  and  control  of  fuel  cell  water 
and  gas  dynamics.  It  is  important  to  note  that  the  fuel  cell  model 
presented  here  is  not  novel  except  in  relating  cell  flooding  to 
performance.  The  unique  contribution  lies  in  applying  this  sim¬ 
ple  isothermal  model  to  well  approximate  the  dynamic  fuel  cell 
response  under  a  range  of  operating  conditions  by  leveraging 
standard  off-the-shelf  sensors  and  actuators. 

This  paper  is  organized  by  first  presenting  the  experimen¬ 
tal  hardware  in  Section  2  followed  by  the  model  of  gas  and 
water  dynamics  in  Section  3.  The  applied  boundary  conditions 
at  the  gas  channels  and  membrane  surface  are  then  presented 
in  Section  4.  The  impact  of  the  liquid  water  on  cell  voltage  is 
modeled  in  Section  5.  The  parameter  identification  methodol¬ 
ogy  is  presented  in  Section  6.  Finally,  the  model  calibration  and 
validation  results  are  shown  in  in  Sections  7  and  8.  A  list  of  the 
model  parameters  is  given  in  Appendix  A. 

2.  Experimental  hardware 

Experimental  results  are  collected  from  a  24-cell  PEMFC 
stack  which  can  deliver  1.4  kW  continuous  power,  capable 
of  peaking  to  2.5  kW.  The  cell  membranes  are  comprised  of 
GORE  ™PRIMEA®  Series  5620  membrane  electrode  assem¬ 
blies  (ME  As).  The  ME  As  utilize  35  p,m  thick  membranes  with 
0.4  mg  cm-2  and  0.6  mg  cm-2  Pt/C  on  the  anode  and  cathode, 
respectively,  with  a  surface  area  of  approximately  300  cm2. 
The  GDL  material,  which  distributes  gas  from  the  flow  fields 
to  the  active  area  of  the  membrane,  consists  of  double- sided, 
hydrophobic,  version  3  ETek™  ELATs®  with  a  thickness  of 
0.43  mm.  The  flow  fields  are  comprised  of  machined  graphite 
plates  with  gas  channels  that  are  approximately  1  mm  wide  and 
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1  mm  deep.  The  flow  pattern  consists  of  semi- serpentine  pas¬ 
sages  on  the  cathode  (30  channels  in  parallel  that  are  16.0  cm  in 
length  with  two  180°  turns)  and  straight  passages  on  the  anode 
(90  channels  in  parallel  that  are  17.1  cm  in  length). 

The  experimental  hardware,  designed  in  collaboration  with 
the  Schatz  Energy  Research  Center  at  Humboldt  State  Uni¬ 
versity,  is  installed  at  the  Fuel  Cell  Control  Laboratory  at  the 
University  of  Michigan.  A  schematic  of  the  major  experimental 
components  along  with  the  measurement  locations  is  depicted  in 
Fig.  2.  A  computer  controlled  system  coordinates  air,  hydrogen, 
cooling,  and  electrical  subsystems  to  operate  the  PEMFC  stack. 
Dry  pure  hydrogen  is  pressure  regulated  at  the  anode  inlet  to  a 
desired  setpoint.  This  pressure  regulation  system  replenishes  the 
hydrogen  consumed  in  the  chemical  reaction.  For  the  majority  of 
the  operational  time,  the  hydrogen  stream  is  dead-ended  with  no 
flow  external  to  the  anode.  Using  a  purge  solenoid  valve,  hydro¬ 
gen  is  momentarily  purged  through  the  anode  to  remove  water 
and  inert  gases.  Humidified  air  (generated  using  a  membrane- 
based  internal  humidifier)  is  mass  flow  controlled  to  a  desired 


stoichiometric  ratio.  Deionized  water  is  circulated  through  the 
system  to  remove  heat  produced  due  to  the  exothermic  chemi¬ 
cal  reaction.  A  fan  is  used  to  thermostatically  control  (on-off) 
the  stack  outlet  coolant  to  a  desired  temperature.  Measurements 
of  the  dry  gas  mass  flow  rates  supplied  to  the  PEMFC  stack  are 
taken  along  with  the  temperature,  pressure  and  relative  humidity 
in  the  inlet  and  outlet  manifolds. 

Due  to  the  lack  of  a  practical  means  to  directly  measure  the 
accumulation  of  liquid  water  within  a  multi-cell  stack,  consec¬ 
utive  anode  purges  and  cathode  surges  (momentarily  increasing 
the  gas  mass  flow  rates)  were  used  to  indicate  the  presence  of 
liquid  water  in  either  the  anode  or  cathode  channels,  as  shown  in 
Fig.  3.  At  approximately  240  s  the  cathode  was  surged,  causing 
an  increase  in  oxygen  partial  pressure  and  cell  voltage.  How¬ 
ever,  this  momentary  voltage  increase  is  not  sustained  following 
the  surge  and  the  general  voltage  decay  due  to  flooding  in  the 
anode  persists.  Following  an  anode  purge,  the  voltage  quickly 
improves  and  then  gradually  decays  until  the  next  anode  purge 
event  is  initiated.  It  is  important  to  note  that  this  gradual  decay  in 


Fig.  2.  Experimental  hardware  employed  and  measurement  locations.  This  figure  is  modified  from  [52]. 
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Fig.  3.  Experimental  data  showing  impact  of  anode  purging  and  cathode  surging 
on  cell  voltages  at  a  constant  nominal  current  density  of  0.3  A  cm-2,  and  an 
operating  temperature  of  T  =  65  °C.  The  first  subplot  shows  the  24  individual 
cell  voltages  in  thin  lines  along  with  the  average  cell  voltage  with  a  thick  line. 
The  second  subplot  shows  the  anode  and  cathode  inlet  total  pressures. 

cell  voltage  could  be  attributed  to  the  accumulation  of  nitrogen 
in  the  anode  which  would  also  be  expelled  during  and  anode 
purge  event.  However,  during  purge  events  a  significant  mass  of 
liquid  water  can  be  visually  detected  leaving  the  anode.  Thus, 
this  work  focuses  on  the  impact  of  anode  flooding  on  cell  voltage 
and  assumes  nitrogen  is  not  the  culprit. 

3.  Modeling  of  gas  and  water  dynamics 

The  model  of  the  reactant  and  water  dynamics  is  presented 
in  the  following  sections,  describing  the  capillary  transport  of 
liquid  water  and  the  diffusion  of  gases  within  the  GDL,  as  well 
as  the  time  varying  boundary  conditions  at  the  membrane  and 
gas  channel  interfaces.  To  approximate  the  spatial  gradients,  the 
gas  diffusion  layer  was  separated  into  discrete  volumes  using 
standard  finite  difference  techniques. 

The  anode  volume  contains  a  mixture  of  hydrogen  and  water 
vapor,  whereas  the  cathode  volume  contains  a  mixture  of  oxy¬ 
gen,  nitrogen,  and  water  vapor.  The  species  concentrations  in 
the  channel  are  calculated  based  on  the  conservation  of  mass 
assuming  the  channel  is  homogeneous,  lumped-parameter,  and 
isothermal.  Under  load,  we  assume  product  water  is  formed  in 
the  vapor  phase. 

This  product  water  vapor,  combined  with  the  water  vapor 
supplied  with  the  cathode  gas  stream,  is  exchanged  between 
the  anode  and  the  cathode  through  the  hydrophilic  membrane. 
The  protons,  liberated  at  the  anode,  transport  water  to  the 
cathode  through  electro-osmotic  drag,  while  back  diffusion 
transfers  vapor  due  to  a  water  vapor  concentration  gradient 
across  the  membrane.  The  net  flux  of  vapor  through  the  mem¬ 
brane  depends  on  the  relative  magnitudes  of  these  transport 


mechanisms.  Although  there  are  many  efforts  to  experimentally 
quantify  back  diffusion  [30-33],  conflicting  results  suggest  an 
empirically  data-driven  identification  of  water  vapor  diffusion 
might  be  a  practical  approach  to  this  elusive  subject.  Constant 
parameters  have  been  used  to  scale  back  diffusion  models  for 
PEMFCs  with  different  membrane  materials  [27,34].  Using  a 
similar  methodology  as  [27],  in  this  paper  the  membrane  water 
transport  algorithm  employs  a  tunable  parameter  to  scale  the 
membrane  water  diffusion  model  in  Ref.  [31]. 

When  the  production  or  transport  of  water  vapor  overcomes 
the  ability  of  the  vapor  to  diffuse  through  the  GDL  to  the  chan¬ 
nel,  the  vapor  supersaturates  and  condenses.  The  condensed 
liquid  water  accumulates  in  the  GDL  until  it  has  surpassed  the 
immobile  saturation  limit  at  which  point  capillary  flow  will 
carry  the  liquid  water  to  an  area  of  lower  capillary  pressure 
(the  GDL-channel  interface).  Liquid  water  in  the  GDL  occu¬ 
pies  the  pore  space,  reducing  the  diffusion  of  the  reactant  gases. 
However,  we  have  found  that  the  reduction  of  the  reactant  con¬ 
centrations  due  to  the  changes  in  the  gas  diffusivity  alone  is 
not  significant  enough  to  degrade  the  voltage  by  the  magni¬ 
tude  experimentally  observed.  Similar  observations  lead  to  the 
consideration  of  the  reactant  diffusion  in  the  catalyst  layer  [35]. 

We  follow  here  a  different  approach  and  instead  of  adding  the 
catalyst  layer  complexity  to  the  model,  we  consider  the  effects  of 
flooding  on  the  area  available  for  diffusion.  The  water  (in  liquid 
and  vapor  phase)  that  wicks  out  of  the  hydrophobic  GDL  to  the 
channel  ultimately  obstructs  the  area  that  reactants  can  diffuse 
through.  This  effect  is  not  easily  modeled  because  the  GDL 
surface  roughness  makes  it  difficult  to  predict  how  much  GDL 
surface  area  is  blocked  by  a  given  volume  of  liquid  water.  Lor  this 
reason,  we  assume  the  liquid  water  at  the  GDL-channel  interface 
forms  a  layer  of  uniform  thickness.  This  water  layer  spreads 
across  the  surface  of  the  GDL  as  the  volume  of  liquid  water 
in  the  channel  increases,  thus  reducing  the  surface  area,  which 
increases  the  calculated  current  density,  in  turn  lowering  the  cell 
voltage  at  a  fixed  total  stack  current.  In  this  model  the  thickness 
of  the  water  layer  is  an  experimentally  tuned  parameter. 

The  estimation  of  the  average  cell  voltage  is  a  function  of 
the  reactant  concentrations  at  the  surface  of  the  membrane,  the 
membrane  water  content,  temperature,  and  the  calculated  cur¬ 
rent  density  based  on  the  reduced  active  area,  which  in  turn 
is  a  function  of  liquid  water  present  in  the  gas  channel.  There 
are  four  experimentally  tunable  voltage  parameters  which  are 
determined  using  linear  least  squares  for  a  given  set  of  mem¬ 
brane  diffusion  and  water  thickness  parameters.  By  comparing 
the  average  measured  cell  voltage  to  the  model  prediction,  these 
parameters  can  be  re-adjusted  to  match  the  rate  of  decay  and 
magnitude  of  the  voltage  degradation.  This  iterative  process 
allows  all  six  tunable  parameters  to  be  identified. 

3.1.  Summary  of  modeling  assumptions 

In  summary,  the  following  general  assumptions  were  made 
in  developing  the  model  presented: 

•  The  volume  of  liquid  water  within  the  GDL  does  not  restrict 

the  volume  occupied  by  the  gases.  The  authors  in  Ref. 
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[36]  indicated  that  the  diffusion  of  gas  through  the  GDL 
occurs  through  a  hydrophobic  macroporous  structure,  where 
as  the  liquid  water  travels  through  the  non-wet  proofed 
pores  (a  microporous  structure),  implying  that  the  pore  vol¬ 
ume  occupied  by  gases  is  fixed.  Examining  the  time  scale 
decomposition  of  the  reactant  and  water  dynamics  [37],  this 
assumption  primarily  influence  the  liquid  water  dynamics  and 
due  to  the  relatively  small  change  in  liquid  water  volume 
between  the  GDL  sections,  has  a  negligible  impact.  However, 
if  different  boundary  conditions  were  applied  which  signifi¬ 
cantly  modified  the  spatial  distribution  of  liquid  water  in  the 
GDL  sections,  this  assumption  should  be  revisited. 

The  internal  cell  structure  (gas  channel,  GDL  and  membrane) 
is  assumed  to  be  isothermal  and  equal  to  the  time  varying 
coolant  outlet  temperature.  However,  the  gas  inlet  temper¬ 
atures  vary  and  are  used  to  calculate  the  water  vapor  mass 
flow  rates  entrained  with  the  supplied  reactants.  Although  it 
is  true  that  a  multi-cell  stack  with  a  large  active  area  will 
undoubtedly  have  thermal  gradients  within  the  cell  structure 
and  impact  water  transport  [38],  this  assumption  is  adequate 
for  estimating  the  temporal  evolution  in  cell  voltage  experi¬ 
mentally  observed  under  both  flooding  and  drying  conditions, 
as  will  be  shown  in  Section  8.  Accounting  for  dynamic  ther¬ 
mal  states  within  the  gas  diffusion  layer  adds  a  significant 
degree  of  model  complexity  which,  while  useful  for  design, 
may  not  be  appropriate  for  control. 

The  gas  channels  are  treated  as  homogeneous  and  lumped 
parameter.  Additionally,  flow  through  the  GDL  is  modeled 
in  one  dimension  which  neglects  the  difference  in  transport 
mechanisms  for  flow  under  the  ribs  versus  under  the  channels. 
Although  models  do  exist  which  characterize  all  these  com¬ 
plex  phenomena,  the  inclusion  of  this  additional  dimension 
has  a  significant  impact  on  the  number  of  internal  states  in  the 
model. 

The  only  mechanism  for  removing  liquid  water  from  the  gas 
channels  is  through  evaporation.  Although  this  is  a  common 
modeling  assumption,  it  could  result  in  an  underestimation 
of  the  total  mass  of  water  (liquid  and  vapor)  removed  from 
the  anode  during  purges.  The  tuned  model  parameters  may 
compensate  for  this  underestimation  but  the  identified  values 
were  physically  reasonable  and  within  ranges  reported  in  lit¬ 
erature  as  discussed  in  Section  6.  It  has  been  shown  [39]  that 
liquid  water  droplet  instability  and  the  resultant  detachment 
from  the  GDL  to  the  gas  channel  can  be  a  significant  liq¬ 
uid  water  removal  mechanism  at  high  current  density  (high 
gas  velocity).  Therefore,  if  this  model  is  to  be  extended  to 
high  current  density  operation,  this  assumption  should  be 
revisited. 

All  gases  behave  ideally.  The  range  of  system  operating  tem¬ 
peratures  and  pressures  permits  the  assumption  of  ideal  gas 
behavior  for  the  gas  constituents  of  interest. 

Hydrogen,  oxygen  and  nitrogen  molecules  do  not  crossover 
through  the  membrane.  Although  these  thin  polymeric  mem¬ 
branes  permit  the  crossover  of  molecules  when  there  is  a 
concentration  gradient  across  the  membrane  [40],  only  the 
water  crossover  at  steady-state  has  been  considered  in  this 
work  for  the  sake  of  model  simplicity. 


•  Due  to  the  relatively  small  gas  flux  within  the  GDL  at  the 

current  density  range  considered,  the  convective  transport  of 

gas  due  to  bulk  flow  was  neglected. 

3.2.  Nomenclature 

This  section  describes  the  nomenclature  used  throughout  this 
paper.  A  list  of  the  parameters  is  provided  in  Appendix  A,  along 
with  values  and  units.  Time  derivatives  are  denoted  as  d()/d t. 
Spatial  derivatives  through  the  GDL  thickness  in  the  membrane 
direction  (y)  are  denoted  as  3(  )/d y.  In  the  presented  model,  all 
equations  have  SI  units  of  Pa,  N,  m,  kg,  s,  and  J  unless  explicitly 
stated. 

The  symbol  a  is  used  for  water  activity,  c  for  molar  concen¬ 
tration  (mol  m(_3)),  (D)  for  effective  diffusivity  (m2  s-1),  Dw 
for  water  vapor  diffusion  coefficient  (m2  s-1),  E  for  the  theoret¬ 
ical  open  circuit  voltage  (V),  i  for  the  nominal  current  density 
(A  cm-2),  zapp  for  the  apparent  current  density  (A  cm-2),  /o  for 
the  exchange  current  density  (A  cm-2),  I  for  the  total  stack  cur¬ 
rent  (A),  Kx i  for  relative  permeability,  n&  for  electroosmotic  drag 
coefficient  (mol  H20/mol  H+),  N  for  molar  flux  (mol  s-1  m-2), 
p  for  pressure  (Pa),  /?sat  for  the  water  vapor  saturation  pressure 
(Pa),  Rev ap  for  the  evaporation  rate  (mol  s-1  m-3),  s  for  the  frac¬ 
tion  of  liquid  water  volume  to  the  total  volume,  S  for  the  reduced 
liquid  water  saturation,  T  for  temperature  (K),  t/ac t  for  the  acti¬ 
vation  voltage  loss  (V),  f/0hmic  for  the  ohmic  voltage  loss  (V), 
Uconc  for  the  concentration  voltage  loss  (V),  v  for  the  measured 
terminal  cell  voltage  (V),  v  for  the  estimated  terminal  cell  volt¬ 
age  (V),  W  for  mass  flow  rate  (kg  s-1),  v  for  the  mass  fraction, 
and  y  for  the  mole  ratio.  Greek  letters  are  used  where  e  is  for  the 
GDL  porosity,  k  for  membrane  water  content  (mol  H2 O/mol 
SO3-),  0  for  relative  humidity  (0-1),  and  co  for  humidity 
ratio. 

The  subscript  amb  is  used  to  represent  ambient  conditions, 
an  for  anode,  c  for  capillary,  ca  for  cathode,  ch  for  channel, 
ct  for  catalyst,  da  for  dry  air,  dg  for  dry  gas,  e  for  electrode 
(an  or  ca),  fc  for  fuel  cell  stack,  H2  for  hydrogen,  in  for  the 
control  volume  inlet  or  input,  j  as  an  index  for  gas  constituents, 
k  as  an  index  for  discretization  (in  time  or  space),  1  for  liquid 
water,  mb  for  membrane,  N2  for  nitrogen,  O2  for  oxygen,  out 
for  the  control  volume  outlet  or  output,  p  for  pore,  rm  for  return 
manifold,  v  for  water  vapor,  and  w  for  water  (gas  and/or  liquid 
phase). 

3.3.  Liquid  water  capillary  transport 

In  hydrophobic  GDL  material,  as  the  GDL  pore  spaces  fill 
with  liquid  water,  the  capillary  pressure  increases,  causing  the 
water  to  flow  to  adjacent  pores  with  less  water.  This  process 
creates  a  flow  of  liquid  water  through  the  GDL,  resulting  in  an 
injection  of  liquid  into  the  channel.  Applying  the  conservation 
of  mass  to  the  GDL  volume,  the  liquid  water  dynamics,  which 
arise  from  capillary  liquid  water  mass  flow,  W\,  and  the  molar 
evaporation  rate,  Rey ap>  can  be  calculated  by 

ds;  _  /  1  \  dWi  _  flevapMy 

d t  \p\sAfcJ  dy  pi 


(1) 
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where  the  mass  of  liquid  water  in  the  GDL  is  expressed  in  terms 
of  liquid  water  saturation,  s9  which  represents  the  fraction  of  the 
liquid  volume  to  the  pore  volume  (s  =  V\/  Vp),  Afc  the  nominal 
fuel  cell  active  area,  p\  the  liquid  water  density,  Mv  the  molecular 
weight  of  water,  and  s  is  the  GDL  porosity. 

The  flow  of  liquid  water  through  the  GDL  is  a  function  of  the 
capillary  pressure  gradient  [41,42]  described  by 


sAfcp\KKr\ 

W\  — - 

/xi 


(2) 


where  pc  is  the  liquid  water  capillary  pressure,  K  the  absolute 
permeability,  /xi  the  viscosity  of  liquid  water,  and  Kv\  =  S3  is  the 
relative  permeability  of  liquid  water.  The  relative  permeability 
function  suggests  more  pathways  for  capillary  flow  are  available 
as  liquid  water  saturation  increases,  and  is  a  function  of  the 
reduced  liquid  water  saturation,  S ,  shown  by 


S  = 


-  for  S[m  <  s  <  1 

1  —  sim 

0  for  0  <  s  <  S[m, 


(3) 


number  of  moles  of  gas  within  the  pore  volume,  Vp,  where 


(6) 


Diffusion  of  hydrogen  and  water  vapor  occurs  in  the  anode 
GDL  and  the  diffusion  of  oxygen  and  water  vapor  occurs  in 
the  cathode  GDL.  As  a  result,  both  the  anode  and  cathode  gas 
diffusion  can  be  modeled  assuming  binary  diffusion.  It  is  impor¬ 
tant  to  note  that  nitrogen  gas  is  present  in  the  cathode.  As  a 
result,  the  nitrogen  concentration  in  the  channel  is  calculated 
and  assumed  to  the  constant  through  the  GDL  since  it  is  not 
involved  in  the  reduction  reaction  at  the  catalyst.  Ternary  dif¬ 
fusion  must  be  assumed  at  both  the  anode  and  the  cathode  if 
nitrogen  cross-over  were  to  be  considered.  The  total  molar  flux 
is  related  to  the  concentration  gradient,  represented  by 


3  Cj 

nj  =  -(dj)^ 


(7) 


where  {Dj)  is  the  effective  diffusivity  of  the  gas  constituents  in 
the  GDL, 


where  S{m  is  the  value  of  the  immobile  saturation  describing  the 
point  at  which  the  liquid  water  path  becomes  discontinuous  and 
interrupts  capillary  flow.  This  capillary  flow  interruption  occurs 
when  s  <  S{m.  The  results  of  capillary  flow  experiments  using 
glass  beads  as  porous  media  show  that  S[m  =  0.1  [41]. 

Capillary  pressure  is  the  surface  tension  of  the  water 
droplet  integrated  over  the  surface  area.  The  Leverett  J-function 
describes  the  relationship  between  capillary  pressure  and  the 
reduced  water  saturation,  S , 

a  cos  0C  9  o 

pc  =  - _ c  [1.4175  -  2.120 S2  +  1.263S3],  (4) 

v K/e  ^ — — - — v -  -  -- 

J(S) 

where  a  is  the  surface  tension  between  water  and  air,  and  0C  is 
the  contact  angle  of  the  water  droplet  [41]. 

Finally,  the  molar  evaporation  rate  is 


(Dj)  =  Dje 


£-0.11 

1-0.11 


0.785 

(1-S)2, 


(8) 


for  two-dimensional  bulk  diffusion  with  flow  perpendicular  to 
the  GDL  carbon  fibers,  where  Dj  is  the  gas  diffusion  coeffi¬ 
cient.  Porosity,  effective  diffusivity  and  liquid  water  saturation 
for  carbon  Toray®  paper  GDL,  are  modeled  from  [41]. 

Finally,  the  general  temporal  derivative  of  gas  concentration 
as  a  function  of  the  local  molar  flux  gradient  and  the  local 
reaction  rate,  Rj,  of  the  particular  gas  species  forms  a  partial 
differential  equation  (PDE), 


dcj_ 

dt 


dy 


+  R  i 


(9) 


where  (7)— (9)  are  combined  to  yield  a  second  order  PDE. 


R 


evap  —  Y 


P  sat  P\ 

RT 


(5) 


where  y  is  the  volumetric  condensation  coefficient  [41],  R  the 
ideal  gas  constant,  T  the  temperature,  pSSLt  the  water  vapor  sat¬ 
uration  pressure  which  itself  is  a  function  of  temperature,  and 
py  is  the  water  vapor  partial  pressure.  When  the  partial  pressure 
of  water  vapor  is  greater  than  the  saturation  pressure,  R&y ap  is 
negative,  representing  the  condensation  of  water.  A  logical  con¬ 
straint  must  be  included  such  that  if  no  liquid  water  is  present 
(s  <  0)  and  the  saturation  pressure  is  greater  than  the  water  vapor 
pressure,  then  water  can  not  be  evaporated  (R&y ap  =  0). 


3.4.  Gas  species  diffusion 

The  diffusion  of  gas  species  in  the  GDL  is  a  function  of  the 
concentration  gradient,  transferring  gas  from  regions  of  higher 
concentration  to  regions  of  lower  concentration.  The  molar  con¬ 
centration  of  gas  species  j  is  denoted  cj  and  is  a  function  of  the 


4.  Boundary  conditions 

The  membrane  and  gas  channels  serve  as  time-varying 
boundary  conditions  for  the  GDL  model.  This  section  presents 
the  application  of  mass  conservation  in  the  channel  as  well  as  the 
model  for  the  water  vapor  exchange  between  the  anode  and  cath¬ 
ode  through  the  membrane.  It  is  important  to  remember  that  the 
spatial  gradients  within  the  GDL  are  approximated  with  finite 
difference  equations.  A  variable  taken  from  a  GDL  section  that 
is  adjacent  to  the  boundary  of  interest  will  be  denoted  by  ^(1) 
or  \I/(L ),  where  (L  =  3)  indicates  the  section  next  to  the  gas 
channel  and  (1)  indicates  the  section  next  to  the  membrane. 

Each  gas  diffusion  layer  is  separated  into  (L  =  3)  discrete 
volumes,  shown  in  Fig.  4,  to  approximate  the  solution  of  (1)  and 
(9)  for  each  of  the  constituents  in  the  GDL.  Spatial  discretiza¬ 
tion  of  the  GDL  yields  eighteen  coupled  ordinary  differential 
equations  (ODEs),  describing  the  gas  constituent  concentrations 
and  liquid  water  saturation,  that  approximate  the  solution  of  the 
original  PDEs. 
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Fig.  4.  Discretization  of  the  gas  diffusion  layers.  The  direction  of  the  assumed  mass  flow  rate  is  indicated  with  a  solid  arrow.  The  dashed  arrow  is  used  to  indicate 
periodic  mass  flow  rates. 


4.1.  Membrane  boundary  conditions 


The  reaction  at  the  catalyst  surface  of  the  membrane  results 
in  a  loss  of  hydrogen  and  oxygen  at  the  anode  and  cathode, 
respectively.  These  fluxes,  are  used  in  the  calculation  of 

the  molar  flux  spatial  gradients,  described  by 


N 7,e(0)  = 


sAfc2^F 


with  §  = 


1  for  j  =  H2 

2  for  j  =  O2 


(10) 


where  I  is  the  total  current  drawn  from  the  stack  and  F  is 
the  Faraday  constant.  The  molar  flux  of  water  vapor  at  the 
GDL-membrane  boundary,  Av?e(0),  is  influenced  by  the  gen¬ 
eration  of  water  vapor  at  the  cathode  membrane  surface  as  well 
as  the  flow  of  water  vapor  through  the  membrane,  such  that 


1 

AVan(0)  —  —  Ft vmb , 
£ 


(Ha) 


Afv,ca(0)  =-(  -1—  +  Av,mb)  .  (lib) 

£  \2FAfc  J 

Note,  a  scaling  factor  of  1  /s  is  used  here  to  ensure  that  the  water 
vapor  mass  flow  rate  through  the  membrane  is  equal  to  the  mass 
flow  rate  entering  the  GDL  at  the  membrane  boundary. 

The  water  content  of  the  membrane  influences  the  membrane 
vapor  transport  which  establishes  a  time-varying  boundary  con¬ 
dition  for  both  the  anode  and  the  cathode.  These  membrane 
properties,  described  in  Ref.  [30],  are  assumed  to  be  invari¬ 
ant  across  the  membrane  surface.  The  spatial  variation  of  water 


vapor  throughout  the  membrane  is  neglected  due  to  the  signif¬ 
icant  difference  in  thickness  between  the  GDL  (432  p,m)  and 
the  membrane  (35  p,m).  It  is  important  to  note  that  the  mem¬ 
brane  transport  properties  presented  in  this  section  are  taken 
from  experimental  work  conducted  at  steady- state.  Non  steady- 
state  phenomena,  such  as  membrane  swelling  and  hysteresis, 
could  be  added  in  the  future  to  improve  model  fidelity. 

As  with  the  other  volumes,  the  membrane  is  considered  to  be 
homogeneous  and  lumped  parameter.  The  flux  of  water  vapor 
through  the  membrane,  Av?mb,  accounts  for  the  effects  of  both 
back-diffusion  and  electro-osmotic  drag,  suggested  by  Springer 
etal.  [30], 


A7  i  ^  (Cv,ca,mb  G,an,mb) 

Ay, mb  —  nd~  —  OiwL)w  , 

F  Gib 


(12) 


where  i  is  the  nominal  fuel  cell  current  density  (7/Afc),  n&  the 
electro-osmotic  drag  coefficient,  Dw  the  membrane  water  vapor 
diffusion  coefficient,  and  Gib  is  the  membrane  thickness.  The 
parameter  aw  is  a  tunable  parameter  that  will  be  identified  using 
experimental  data.  The  convective  water  transport  mechanisms 
suggested  in  Refs.  [7,38,43]  are  neglected  due  to  the  relatively 
small  water  pressure  gradients  at  these  operating  conditions. 

The  electro-osmotic  drag  coefficient,  described  by  Springer 
et  al.  [30],  is  calculated  using 


= 


2.5A.mb 


(13) 
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where  the  membrane  water  content,  A,mb  is  defined  as  the  ratio 
of  water  molecules  to  the  number  of  charge  sites. 

The  water  vapor  concentration  in  the  electrode  at  the  mem¬ 
brane  surface  is 

c«„b  =  J,  (14) 

47mb,dry 


where  pmb,dry  is  the  membrane  dry  density,  Mmb,dry  is  the  mem¬ 
brane  dry  equivalent  weight,  and  ke  is  the  membrane  water 
content  at  the  surface  of  the  membrane  next  to  either  the  anode 
or  cathode  GDL. 

The  water  vapor  diffusion  coefficient  for  a  perflourinated 
ionomeric  membrane,  Nation®  117,  was  determined  at  25  °  C 
by  Fuller  [3 1]  by  applying  a  mass  balance  to  determine  the  water 
vapor  flux  through  the  membrane,  resulting  in 


Dw  =  3.5  x  10“6 


-2436 

T 


(15) 


Two  different  cubic  polynomials  were  presented  by  Springer 
et  al.  [30]  and  Hinatsu  et  al.  [44]  to  relate  water  activity  to 
membrane  water  content  at  30  °  C  and  80  °C,  shown  as 


kf°c  =  0.043  +  17.81a/  -  39.85a]  +  36.0a]  (16a) 

kf°c  =  0.300  +  10.8a;  -  16.0a]  +  14.1a]  (16b) 


where  a  is  the  water  activity  and  the  subscript  j  is  used  here 
to  distinguish  between  the  anode  or  cathode  membrane  surface 
and  within  the  membrane  itself,  j  e  {an,ca,mb}.  To  estimate  the 
water  content  at  intermediate  temperatures  and  sub-saturated 
conditions  [45],  suggested  a  linear  interpolation  between  the 
two  uptake  isotherms  shown  in  (16),  such  that 


kj  =  (k 


80  °C 


■kf°C) 


T  -  303 
353  -  303 


30  °C 


(17) 


It  is  important  to  note  that  these  two  uptake  isotherms  are 
applicable  only  when  water  is  in  the  vapor  phase  (a j  <  1). 

In  Ref.  [30]  it  was  shown  that  a  membrane  equilibrated  with 
liquid  water  has  a  water  content  of  A,  =  16.8  at  80  °C,  which  dif¬ 
fers  from  the  water  content  when  the  membrane  is  equilibrated 
with  a  saturated  vapor.  It  was  further  indicated  that  the  water 
content  is  sensitive  to  temperature  when  equilibrated  with  liquid 
water,  but  assumed  to  be  a  linear  relationship  between  [k  =  14, 
a  —  1]  and  \k  =  16.8,  a  —  3]  regardless  of  temperature,  due  to 
a  lack  of  data  regarding  the  membrane  equilibration  for  water  in 
both  the  liquid  and  vapor  phase.  Similarly,  we  assume  a  linear 
relationship  between  the  membrane  water  content  when  equili¬ 
brated  with  water  vapor,  shown  in  (17)  for  aj  —  1,  and  the  value 
of  k  —  16.8  at  a  —  3  published  by  Springer  et  al.  [30],  such  that 


16.8  -  A/;=l  \ 


(aj  —  V)  +  kaj  1 


(18) 


for  1  <  aj  <  3.  Further  experimental  results  from  [46]  and  [44] 
provided  data  regarding  the  temperature  sensitivity  of  the  mem¬ 
brane  water  content  equilibrated  with  liquid  water.  However, 
fitting  this  data  points  to  a  non-monotonic  behavior  of  k  =  f(a ), 
at  some  temperatures  within  the  operating  range  of  the  PEMFC, 


during  the  transition  between  water  in  the  vapor  and  liquid 
phases  (a  =  1),  hence  this  relationship  is  not  considered  in  this 
work. 

Finally,  the  membrane  water  activity  is  assumed  to  be  the 
average  between  the  anode  and  cathode  water  activities  (defined 
by  the  GDL  sections  closest  to  the  membrane  surface),  described 
by 


^mb  — 


^an(l)  +  ^ca(l) 
2 


and 


fled)  = 


Pv,e(  1) 
Psat 


(19) 


where  pv,e( 1 )  is  the  water  vapor  pressure  in  the  GDL  layer  next  to 
the  membrane,  calculated  using  the  water  vapor  concentrations. 

Note  here,  it  is  assumed  that  reactant  molecules  do  not  trans¬ 
fer  through  the  membrane  between  the  anode  and  the  cathode. 
Additionally,  only  water  vapor  can  penetrate  the  membrane,  not 
liquid  water,  implying  lTi,e(0)  =  0. 

In  summary,  the  water  vapor  partial  pressures  in  the  GDL 
section  closest  to  the  membrane  surfaces  are  used  to  determine 
the  water  activity  in  the  first  GDL  section,  which  is  assumed  to 
be  equal  to  the  membrane  water  activity  at  the  membrane-GDL 
interface.  These  two  membrane  water  activities  are  averaged 
to  calculate  the  lumped  membrane  water  activity,  which  influ¬ 
ence  diffusion  and  electro-osmotic  drag.  Finally,  the  net  water 
vapor  flux  is  calculated,  given  diffusion,  drag  and  the  water  vapor 
concentrations  at  the  membrane  surfaces. 


4.2.  Boundary  conditions  at  the  cathode  channel 

The  concentration  of  oxygen  and  water  vapor  in  the  cath¬ 
ode  channels,  co2,ca(f-  +  1)  and  cV;Ca(L  +  1),  are  used  for  the 
calculation  of  the  gas  concentration  gradient  for  the  GDL  sec¬ 
tion  next  to  the  channels,  (3c;?ca/3y)(L).  Mass  conservation 
for  the  gas  species  in  the  cathode  is  applied  using  the  cathode 
inlet  conditions  as  inputs,  requiring  measurements  of  the  dry  air 
mass  flow  rate,  TTda,Ca,in,  temperature,  rca? in,  total  gas  pressure, 
Pea, in,  and  humidity,  0Ca,in,  along  with  the  cathode  outlet  pres¬ 
sure,  Pea, out-  After  completing  several  experiments  under  a  range 
mass  flow  rates  and  temperatures,  it  was  found  that  the  cathode 
inlet  total  gas  flow  was  fully  humidified  and  the  cathode  outlet 
total  pressure  was  approximately  atmospheric,  motivating  the 
assumptions  that  0Ca,in  =  1  and  pca,out  =  Patm- 

The  mass  flow  rate  of  the  individual  gas  species  supplied  to 
the  cathode  channel  are  calculated  as  follows: 


W02  ,ca,in  —  -^02,ca,infkda,ca,in,  WN2 

,ca,in  —  -^N2,ca,inlkda,ca,in, 

fkv,ca,in  =  ^ca,in  fkda,ca,in  >  (70) 


where  the  humidity  ratio,  co,  is  generally  defined  by 


My  0Psat 
^dg  P  -  0Psat 


(21) 


for  a  gas-water  vapor  mixture,  with  the  mass  fraction  of  oxygen 
and  nitrogen  in  the  dry  air  (da)  defined  as  xo2  =  yo2  Mo2  /  M da 
and  xn2  =  (1  -  yo2)^N2/^da,  where  M&d  =  yo2Mo2  +  (1  - 
yo2)^N2  and  yo2  is  the  oxygen  mole  fraction  in  dry  air. 
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The  gas  species  mass  in  the  cathode  channel  are  balanced  by 
applying  mass  continuity: 


dmo2;Ca(£  +  1) 
d  t 


—  Wo2, ca, in  Wo2,ca,out  H”  Wo2, ca(7>), 

dmN2,ca(L  +  1) 


d  t 

d^w,ca(^  + 1) 
d  t 


-=Wn. 


2,ca,m 


WN2,ca,out, 


—  Wv?ca,  in  Wv^outH- Ww?ca(L). 

(22) 


The  cathode  channel  pressure  is  calculated  applying  Dalton’s 
law  such  that 


rr  \  i\  _  ^  (  m02,ca(L  +  1)  (  m^2^ca(L  +  1) 
Pca\P  “r  -U  —  77  I  —  I  — 

Tca  V  Mq2  MN2 


+min 


Psaii 


Pd  g,ca(^  +  l) 

RTmw,ca(L  +  1) 

MyVca 


(23) 


Pv,ca(^+l) 


Although  in  the  physical  system  the  cathode  air  mass  flow 
rate  may  be  responsible  for  removing  some  liquid  water  from 
the  cathode  channel,  for  modeling  purposes  it  is  assumed  that 
all  water  exiting  the  cathode  is  in  the  form  of  vapor. 

The  mass  flow  rate  of  gases  exiting  the  cathode  are  calculated 
as 


the  liquid  water  mass  flow  rate  between  the  GDL-channel  inter¬ 
face,  WijCa(L).  Within  the  channel,  the  volume  of  liquid  water 
is  assumed  to  be  negligible  compared  with  the  total  channel 
volume,  motivating  this  assumption  that  Scsl(L  +  1)  =  0. 


4.3.  Boundary  conditions  at  the  anode  channel 


Similarly  to  the  cathode,  the  inputs  for  the  anode  calculations 
are  the  measured  anode  inlet  conditions  including  the  dry  hydro¬ 
gen  mass  flow  rate,  WH2,an,in,  the  supply  manifold  temperature, 
Tan, in,  the  total  pressure,  pa n,in,  and  the  relative  humidity,  0an,in- 
Dry  hydrogen  is  supplied  to  the  anode,  as  a  result  0an,in  =  0. 
The  resulting  mass  balances  for  hydrogen  and  water  are 


d^H2,an(T>  +  1) 
d  t 

draW,an(T  +  1) 
d  t 


—  WH2,an,in  WH2,an,out  f^H2,an 
=  Wv,an,in  —  Wv,an,out  —  Ww5an(L). 


(L\ 


(26) 


The  dry  hydrogen  inlet  mass  flow  rate,  WH2,an,in  = 
£an,in(/?an,in  —  Ptm(L  +  1)),  is  controlled  with  a  pressure  regu¬ 
lator  to  maintain  a  constant  anode  inlet  total  pressure.  Because 
the  hydrogen  supplied  to  the  anode  is  dry,  the  vapor  mass  flow 
rate  is  assumed  to  be  zero  ( Wv,an,in  =  0).  In  calculating  the  anode 
total  channel  pressure,  both  the  partial  pressures  of  hydrogen  and 
water  vapor  must  be  estimated  such  that, 


Wca,out  —  kca(Pca(L  +  1)  Pea, out), 

1 

f^da,ca,out  —  z  .  Wca?0ut, 

I  +  <^ca,out 

Wo2  ,ca,out  —  -^02,ca,ch  Wda,ca,out, 
Wv,ca,out  =  Wca?0ut  —  Wa;Ca;0ut, 
WN2,ca,out  =  (1  —  -^02,ca)f^da,ca,out, 


(24) 


RT 

Pan(L  +  1)  —  —  —  ^2  H2  (L  +  1) 

Mu2  Van 


+min 


/?H2,an(£  +  l) 

RTm 


Psat, 


:w,an 


L  (L  +  1) 


MvVaT 


(27) 


where  kc&  is  an  orifice  constant  found  experimentally.  Although 
the  mole  fraction  of  oxygen  at  the  cathode  inlet  is  assumed  to 
be  constant,  yo2,ca,in  =  0.21,  the  mole  fraction  of  oxygen  in  the 
channel  (driving  the  outlet  mass  flow  rates)  is  dependent  upon  the 
oxygen  mass  (pressure)  state  in  the  channel,  such  that  yo2,ca  = 
P02,cJ  Pea- 

Finally,  the  oxygen  and  total  water  mass  flow  rates  between 
the  GDL  and  the  channel,  Wo2,ca(£)  and  WW;Ca(L),  must  be  cal¬ 
culated  to  solve  the  mass  conservation  equations  shown  in  (22). 
The  oxygen  mass  flow  through  the  GDL-channel  interface  is  a 
function  of  the  oxygen  molar  flux,  No2  (L).  The  total  water  mass 
flow  rate,  Ww?ca(L),  exchanged  between  the  GDL  and  channel 
is  a  function  of  the  liquid  water  mass  flow,  Wi?ca(L),  and  the 
water  vapor  flux,  Av,Ca-  Both  the  oxygen  and  total  water  mass 
flow  rates  are  described  by 

Wo2,ca(£)  =  No2(L)Mo2£Afcnce\\s, 

WW;Ca(L)  =  (Wi>ca(L)  +  Ny^(L)My£Afc)  nCQ ns,  (25) 

where  the  assumption  5’ca(L  +  1)  =  0  is  employed  in  the  cal¬ 
culation  of  the  reduced  water  saturation  gradient  to  determine 


The  total  mass  flow  rate  leaving  the  anode  channel,  Wan, out, 
exists  only  during  an  anode  gas  purge  to  remove  both  water,  and 
unfortunately,  hydrogen.  The  equations  quantifying  the  hydro¬ 
gen  and  water  vapor  mass  flow  rates  leaving  the  anode  channel 
are  expressed  as 

Wan,out  =  ^an,out(Pan(T  +  1)  —  pan,out), 

1 

WH2,an,out  —  “  '  Wan?0ut, 

1  +  &>an,out 

Wv,an,out  =  Wan?out  —  Wn2,an,out-  (28) 

Similarly  to  the  cathode,  the  gas  and  liquid  water  mass  flow 
rates  between  the  GDL  and  channel  are  calculated  by 

WH2,an(T)  =  Nn2(L)Mn2£AfcnCQ\\s, 

Ww?an(L)  =  (Wi}an(L)  +  Afy5an  (L)MV£  A  fc)ftceiis  (29) 

where  the  assumption  San(L  +  1)  =  0  is  employed  in  the  calcu¬ 
lation  of  the  reduced  water  saturation  gradient  to  determine  the 
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liquid  water  mass  flow  rate  between  the  GDL-channel  interface, 
Wi,an(£). 

The  calculation  of  the  mass  flow  rates  leaving  the  anode 
channel  depends  on  the  measurement  of  the  anode  outlet  total 
pressure,  pan,ouu  shown  in  (28).  The  anode  outlet  pressure 
can  also  be  estimated  using  a  similar  approach  as  presented 
for  the  anode  channel  and  documented  in  Ref.  [26],  where 
Wan,rm  =  &an,rm(Pan,out  -  Pamb),  resulting  in  the  addition  of  two 
states  (hydrogen  and  water  mass  in  the  return  manifold). 


Section  7,  this  tuned  parameter  is  similar  to  that  experimentally 
determined  in  Ref.  [6]. 

Once  the  apparent  current  density  is  calculated  it  is  used, 
together  with  the  partial  pressure  of  the  reactants  in  the  anode 
and  cathode  GDL  sections  next  to  the  membrane,  to  determine 
the  average  cell  voltage.  The  average  cell  voltage,  v,  is  equal 
to  the  theoretical  open  circuit  voltage,  E ,  minus  the  activation, 
T/ act?  and  ohmic,  £/0hmic>  losses  such  that 

n  =  E  —  C/act  —  f^ohmic-  (33) 


5.  Output  voltage  equation 


In  this  section,  the  voltage  equation  is  presented  as  a  map¬ 
ping  from  the  apparent  current  density,  reactant  concentrations, 
temperature  and  membrane  humidity  conditions.  All  units  for 
current  density  used  throughout  the  presentation  of  the  voltage 
model  are  given  in  A  cm-2  for  consistency  with  other  published 
models. 

Once  anode  flooding  occurs,  we  associate  the  resulting  volt¬ 
age  degradation  with  the  accumulation  of  liquid  water  mass  in 
the  anode  channel, 


mi>an(L  +  1)  =  max 


0,  Wlw, anC^  T  1) 


Psat-^v  Van 
RT 


(30) 


where  the  mass  of  water  in  the  anode  channel,  mw?an(L  +  1),  is 
taken  from  (26).  The  accumulated  liquid  water  mass  is  assumed 
to  form  a  thin  film  of  thickness,  tw i,  blocking  part  of  the  active 
fuel  cell  area,  Afc,  and  consequently  increasing  the  apparent 
current  density  [20] 


*app(A  cm  2) 


/(A) 

10, 000Aapp(m2)  ’ 


(31) 


where  the  apparent  fuel  cell  area  Aapp  is  approximated  as 


We  have  assumed  that  the  concentration  voltage  loss  due  to  a 
mass  transport  limitation  at  high  current  density  is  negligible 
as  a  result  of  our  operation  at  relatively  low  current  densities 
(/  <  0.4  A  cm-2). 

The  theoretical  open  circuit  voltage,  if  the  chemical  reaction 
was  a  reversible  process,  varies  with  respect  to  reactant  partial 
pressures  and  temperature  according  to  the  change  in  Gibbs  free 
energy  and  the  Nernst  equation  [50], 

(AH  TAS\  RT  ( PH2,an(l)\/P02,ca(l)\ 

V  2 F  2F  )  +  2F  n{  fa,)1-5  J 

(34) 

where  AS  and  AH  are  the  differences  in  entropy  and  enthalpy 
from  standard  state  conditions,  po  the  standard  pressure,  and  the 
oxygen  and  hydrogen  partial  pressures,  /?o2,ca(l)  and  /?H2,an(l), 
are  located  in  GDL  Section  1  next  to  the  membrane. 

The  activation  overvoltage  accounts  for  the  energy  required 
to  drive  the  chemical  reaction  (a  deviation  from  equilibrium),  as 
well  as  the  loss  current  density  resulting  from  the  transport  of 
molecular  hydrogen  from  the  anode  to  the  cathode  through  the 
membrane.  The  total  activation  voltage  loss  was  parameterized 
according  to  [51],  such  that 


aPP 


—  -Afc 


2mun(L  +  1) 
n  cells  Pi  flvl 


(32) 


RT 

Aact  =  A  i  — —  In 
E 


Gpp  +  floss 

io 


(35) 


The  scaling  factor  of  2  in  (32)  was  used  to  account  for  the  fact  that 
one  half  of  the  surface  area  at  the  GDL-channel  interface  is  occu¬ 
pied  by  channel  ribs,  which  reduces  the  area  available  for  the 
formation  of  a  liquid  water  film.  This  methodology  for  relating 
the  accumulation  of  the  liquid  water  in  the  channel  to  a  restricted 
active  area  was  first  proposed  in  Ref.  [20]  and  a  similar  method¬ 
ology  was  employed  by  Hernandez  et  al.  [29].  Some  models  that 
deal  with  cathode  flooding,  however,  propose  an  increased  cur¬ 
rent  density  due  to  the  water  accumulation  in  the  catalyst  layer  at 
the  GDL-membrane  interface  [47].  Ongoing  experimental  work 
from  many  researchers  has  focused  on  quantifying  this  accumu¬ 
lation  of  liquid  water  using  direct  visualization  [48]  or  neutron 
imaging  techniques  [3,6,49]. 

The  thickness  of  this  water  layer,  tw\  is  a  tunable  parameter 
that  impacts  the  rate  at  which  the  active  area  is  reduced  and  in 
turn  the  rate  of  voltage  decay  as  the  liquid  water  accumulates. 
Note  that  the  notion  of  apparent  current  density,  influenced  by  tw\ 
in  the  gas  channel,  is  a  simplification  of  the  flooding  phenomena 
that  nevertheless  captures  the  experimentally  observed  dynamic 
voltage  behavior  of  a  multi-cell  stack  under  a  range  of  condi¬ 
tions  including  both  flooding  and  non-flooding.  As  shown  in 


where  K\  is  a  tunable  parameter  representing  the  reciprocal  of 
the  charge  transfer  coefficient,  /i0Ss  the  loss  current  density  due 
to  hydrogen  crossover,  zapp  the  apparent  current  density  that  is 
a  function  of  the  reduced  active  area  due  to  the  accumulation  of 
liquid  water  at  the  GDL-channel  interface  from  (31),  and  io  is 
the  exchange  current  density  which  is  a  function  of  the  reactant 
partial  pressure  and  temperature  [51],  expressed  as 


(36) 


where  AT  2  and  A3  are  tunable  parameters,  Ec  the  activation 
energy  for  oxygen  reduction  on  Pt,  and  To  is  the  reference 
temperature. 

The  ohmic  voltage  loss  is  dominated  by  the  membrane  con¬ 
ductivity  as  well  as  the  contact  and  bulk  electrical  resistance  of 
the  conductive  materials.  This  loss  was  shown  experimentally 
in  Ref.  [30]  to  have  the  following  functional  form, 


Lohmic  —  A 4 


flnb 


-1268 


(flllAmb  ^12) 


fJ-  -  L\ 

y 303  T J 


*app> 


(37) 
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where  K4  is  a  tunable  parameter,  tm b  is  the  membrane  thickness, 
bn  and  b\2  the  experimentally  identified  parameters  from  [30], 
and  Amb  is  the  membrane  water  content  from  (16)— (1 8). 


-K4 


hnb 


-1268 


(^ll^mb  -  &12) 


(J L 

(  303 


O'app  4"  Ooss)-  (40) 


6.  Parameter  identification  approach 


Lacking  a  practical  experimental  means  to  measure  the  spa¬ 
tial  distribution  of  water  mass  in  the  anode  and  cathode  of  a 
large  multi-cell  stack  for  the  use  of  online  control,  the  lumped- 
parameter  two-phase  flow  model  developed  here  is  indirectly 
calibrated  and  validated  through  model  prediction  of  the  effects 
of  flooding  on  cell  voltage.  A  reasonably  wide  variation  in  the 
experimental  operating  conditions  have  been  examined,  includ¬ 
ing  both  flooding  and  non-flooding  conditions,  to  ensure  that 
the  model  adequately  estimates  the  relationship  between  GDL 
flooding  and  cell  voltage  degradation.  The  range  of  operating 
conditions  examined  is  limited  due  to  our  operation  with  a  stack, 
not  a  single  cell,  and  our  desire  to  minimize  cell  to  cell  voltage 
variations  [52]. 

There  exists  two  sets  of  model  parameters  which  must  be 
either  calibrated  or  tuned.  The  calibrated  parameters  are  based 
on  the  fuel  cell  hardware  specifications  and  are  listed  in  Table  1 
with  values  provided  in  Appendix  A.  These  parameters  may 
require  additional  experiments  to  determine,  such  as  the  orifice 
constants  describing  the  back  pressure  flow  characteristics  for 
each  gas  channel. 

The  two  water  related  tunable  parameters  that  require  exper¬ 
imental  identification  are  the:  scaled  “stack-level”  membrane 
back  diffusion,  aw,  of  (12),  and  thickness  of  liquid  water  layer 
accumulating  at  the  GDL-channel  interface,  tw\,  of  (32).  Addi¬ 
tionally,  there  are  four  tunable  parameters  K\-  K 4  associated 
with  the  output  voltage  in  (33)-(37).  Although  the  water  related 
parameters  do  not  appear  linearly,  the  voltage  equation  can  be 
rearranged  such  that  each  of  the  tunable  K’  s  is  linear  in  the 
coefficient, 

v=E-Kl—  f  ln(iapp  +  /loss)  +  -j  (  =  -  jr  J  J  (38) 


RT  RT  ( p(h  ca(l) 

+  ln(tf2)  *i  —  +  £3  Ki  —  In  TO2'caV 
F  F  \  po 


(39) 


Table  1 


Parameters  required  based  on  PEMFC  stack  specifications 


Afc 

Fuel  cell  nominal  active  area 

K 

Absolute  permeability 

Aim, dry 

Membrane  dry  equivalent  weight 

n 

Number  of  cells  in  the  stack 

%dl 

GDL  thickness 

finb 

Membrane  thickness 

tan 

Total  anode  channel  volume 

Vca 

Total  cathode  channel  volume 

6 

GDL  porosity 

Pm, dry 

Membrane  dry  density 

^an 

Anode  orifice  constants 

kca 

Cathode  orifice  constants 

Note:  values  for  these  parameters  are  provided  in  Appendix  A. 


Given  a  set  of  values  for  aw  and  twb  the  voltage  parameters 
were  identified  using  linear  least  squares  to  minimize  the  dif¬ 
ference  between  the  measured  average  cell  voltage,  v,  and  the 
modeled  cell  voltage,  v , 

/^exP 

[S(t)  -  6(t)]t[i)(t)  -  t>(r)]  dr,  (41) 

over  the  experimental  testing  time,  feXp-  The  statistics  associated 
with  the  estimation  error  were  examined  over  a  range  of  [aw ,  fwi] 
pairs  to  find  the  locally  optimal  [aw,  twi ]  combination  and  the 
resulting  K  values. 

Note  here,  there  are  24  individual  cell  voltages  being  mea¬ 
sured.  The  average  and  median  cell  voltages  exhibit  similar 
dynamics  with  a  relatively  small  difference  in  voltage  between 
them.  However,  there  is  a  significant  difference  in  the  magnitude 
and  deviation  between  the  minimum  and  maximum  cell  volt¬ 
ages.  As  a  result,  the  use  of  either  the  minimum  or  maximum 
cell  voltages  for  parameter  tuning  results  in  an  underestimation 
or  overestimation  of  the  degree  of  flooding.  For  these  reasons, 
the  average  cell  voltage  is  used  for  model  tuning. 

7.  Model  calibration 

Experimental  calibration  data  were  collected  for  a  range  of 
nominal  stack  current  densities  from  /  =  0  to  300  mA  cm-2,  air 
stoichiometries  of  250%  and  300%,  and  coolant  outlet  tempera¬ 
tures  from  45  to  63  °C,  at  an  anode  inlet  total  pressure  of  1 .2  bar, 
as  shown  in  Fig.  5.  A  polarization  curve  ( I-V)  was  conducted 
at  approximately  70  min,  at  which  time  the  purge  events  were 
temporarily  disabled. 

The  purge  events  were  scheduled  to  occur  every  180  s  for  a 
duration  of  1  s.  During  purge  events,  the  purge  solenoid  valve 
was  momentarily  opened,  exposing  the  anode  outlet  manifold 
to  ambient  pressure.  As  a  result  of  this  decreased  anode  total 
pressure,  the  manual  pressure  regulator,  which  tries  to  maintain 
its  downstream  pressure,  increased  the  hydrogen  mass  flow  rate 
through  the  system.  Following  the  closure  of  the  purge  solenoid 
valve,  small  spikes  in  pressure  occur  as  the  pressure  regulator 
readjusted  its  delivery  pressure. 

As  shown  in  Fig.  5,  the  initial  coolant  outlet  temperature 
setpoint  was  50°  C  and  then  changed  to  60°  C  at  approxi¬ 
mately  185  min.  Thermostatic  controllers  were  used  to  control 
the  heat  exchanger  fans  to  regulate  the  coolant  outlet  tempera¬ 
ture.  As  these  fans  were  cycled,  oscillations  in  temperature  were 
induced. 

The  standard  deviation  in  the  cell  voltage  measurements  was 
greater  at  high  current  density  (300  mA  cm-2)  than  at  low  cur¬ 
rent  density  [52].  This  increased  uncertainty  at  high  current 
density,  seen  in  Fig.  5,  was  due  to  both  the  increased  differ¬ 
ence  in  the  cell  to  cell  voltage  variation  as  well  as  the  increased 
excursions  in  cell  voltage  between  anode  purges.  Moreover,  at 
high  current  density  the  cell  with  the  minimum  voltage  exhib¬ 
ited  greater  voltage  excursions  between  anode  purges  than  the 
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Fig.  5.  Experimental  measurements  used  for  model  calibration.  The  first  subplot 
shows  the  24  individual  cell  voltages  along  with  black  dots  at  600  mV  which 
illustrate  the  portions  of  the  data  set  used  for  calibration.  The  second  subplot 
shows  the  nominal  current  density.  The  third  subplot  is  the  anode  and  cathode 
inlet  total  pressures.  The  fourth  subplot  is  the  temperature  of  the  water  coolant 
leaving  the  stack. 

cell  with  the  maximum  voltage.  However,  the  mean  and  median 
voltages  had  similar  dynamic  and  steady- state  responses. 

For  the  purposes  of  model  calibration,  a  portion  of  the  cali¬ 
bration  data  set  was  selected  to  include  a  range  of  both  transient 
and  “steady-state”  operating  conditions.  This  portion  of  data  is 
indicated  with  a  black  v  in  the  voltage  plot  shown  in  Fig.  5.  Data 
at  open-circuit  were  avoided  due  to  the  high  uncertainty  associ¬ 
ated  with  operation  at  open-circuit  voltage  [52].  The  identified 
parameters  resulting  in  the  smallest  mean,  maximum  and  stan¬ 
dard  deviation  in  the  estimation  error  over  the  set  of  aw  £  [7,  12] 
and  tw\  e  [0.12mm,0.14mm],  while  still  capturing  the  voltage 
response  during  flooding  conditions,  are  shown  in  Table  2. 

The  first  tunable  voltage  parameter,  K i,  which  scales  the 
total  activation  overvoltage  has  a  local  minimum,  whose  value  is 
dependent  on  tw\,  near  aw  =  10.5.  The  second  and  third  tunable 
voltage  parameters,  K 2  and  Ky,  influence  the  exchange  current 
density  and  tend  to  increase  as  increases  or  tw  decreases. 


Table  2 


Experimentally  identified  parameter  values 


Parameter 

Tuned  value 

Ki 

1.00 

K2  (pA) 

1.24 

k3 

2.05 

ka 

3.40 

Q!w 

10.0 

twl  (mm) 

0.14 

The  fourth  tunable  voltage  parameter,  K 4,  scales  the  ohmic 
overvoltage  and  decreases  as  aw  increases  or  tw  decreases. 

Using  the  identified  parameters,  the  model  was  simulated  to 
produce  voltage  estimations  for  the  entire  calibration  data  set. 
Fig.  6  shows  the  model  estimation  at  300  mAcm-2  between 
180  and  200  min.  The  second  subplot  compares  the  nominal 
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Fig.  6.  Model  calibration  results.  The  first  subplot  shows  the  24  individual  cell 
voltages  in  thin  faint  colored  lines  with  the  average  cell  voltage  in  a  thick  solid 
blue  line  and  the  model  estimated  average  cell  voltage  in  a  thick  dotted  red  line. 
The  second  subplot  shows  the  nominal  and  apparent  current  densities.  The  third 
subplot  is  the  anode  and  cathode  inlet  total  pressures.  The  fourth  subplot  is  the 
temperature  of  the  water  coolant  leaving  the  stack. 
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current  density,  i  =  7/Afc,  to  the  apparent  current  density,  /app 
from  (31),  based  on  the  apparent  surface  area  that  is  not  blocked 
by  the  liquid  water  film  at  the  GDL-channel  interface. 

As  liquid  water  accumulated  in  the  anode  gas  channels,  the 
apparent  area  decreased,  causing  an  increase  in  the  apparent  cur¬ 
rent  density.  Following  a  purge,  the  liquid  water  was  removed 
and  the  apparent  current  density  returned  to  the  nominal  value. 
Following  some  purges,  not  all  of  the  water  was  removed  from 
the  gas  channels,  causing  the  apparent  current  density  to  remain 
greater  than  the  nominal  current  density.  Since  the  apparent  cur¬ 
rent  density  was  used  to  calculate  the  cell  voltage,  the  estimation 
of  cell  voltage  is  then  sensitive  to  the  degree  of  flooding  in 
the  anode  gas  channels  and  GDL.  The  values  for  the  identified 
parameters,  aw  and  tw\  influence  the  rate  at  which  liquid  water 
accumulates  in  the  gas  channels  (impacting  the  rate  of  decay  in 
voltage  between  purges)  as  well  as  how  much  liquid  water  mass 
accumulates  in  the  gas  channel  (how  much  the  voltage  recovers 
following  a  purge).  When  all  of  the  liquid  water  was  removed 
from  the  gas  channels,  the  cell  voltage  returned  to  approximately 
the  same  value  following  each  purge  event. 

Although  the  voltage  prediction  is  an  indirect  means  for  eval¬ 
uating  the  overall  predictive  ability  of  our  model,  voltage  is  a 
stack  variable  that  combines  the  internal  states  of  the  stack  and 
provides  an  accessible,  cheap,  fast  and  accurate  measurement. 
The  model  accurately  captured  the  trend  of  the  voltage  decay 
and  subsequent  recovery  after  an  anode  purging  event.  Note 
here,  for  the  entire  calibration  data  set,  the  average  estimation 
error  was  2.9  mV,  the  maximum  estimation  error  was  42  mV  and 
the  standard  deviation  in  the  estimation  error  was  3.6  mV. 

In  addition  to  adequately  capturing  the  temporal  evolution 
in  voltage  during  flooding,  the  model  accurately  estimated  the 
reactant  dynamics  during  load  changes.  The  overshoot  in  cell 
voltage  during  a  step  change  up  in  current  from  0.25  A  cm-2 
to  0.3  A  cm-2  at  approximately  183.4  min,  is  shown  in  detail 
in  Fig.  7,  along  with  subsequent  purging  events  near  183.7 
min.  A  decrease  in  the  partial  pressure  of  oxygen  at  the  cath¬ 
ode  membrane  surface  occurs  due  to  volume  filling  dynamics; 
however,  there  was  very  little  deviation  in  the  hydrogen  par¬ 
tial  pressure  during  the  load  change.  As  a  result,  the  reactant 
starvation  occurred  predominantly  on  the  cathode  and  not  the 
anode  under  these  operating  conditions.  Referring  back  to  Fig.  6, 
the  overshoot  in  cell  voltage  at  approximately  198  min  for  a 
step  change  down  in  current  from  0.3  Acm_2to  0.25  Acm-2is 
also  well  approximated.  Note  here,  for  simulations  of  reactant 
dynamics  during  a  load  change  reported  in  Refs.  [23,53],  the 
model  predictions  were  not  compared  with  experimental  data. 

Fig.  8  displays  the  predicted  water  vapor  partial  pressures, 
the  liquid  water  saturation,  and  the  mass  of  liquid  water  accu¬ 
mulating  in  the  anode  channel  during  the  same  load  change  and 
subsequent  purging  events  as  described  previously  for  Fig.  7. 
The  slow  rise  in  the  water  vapor  partial  pressure  was  due  to  the 
increase  in  the  cell  operating  temperature. 

Immediately  following  the  purge  valve  opening,  the  mass 
of  liquid  water  in  the  anode  channel  was  evaporated  into  the 
bulk  gas  stream  (due  to  the  increased  hydrogen  mass  flow  rate 
during  the  purge).  The  volumetric  condensation  coefficient,  y 
in  (5),  influenced  the  non-instantaneous  rate  of  evaporation  of 


Time  (min) 

Fig.  7.  Reactant  dynamics  during  a  load  change  and  purging  event.  The  first  row 
of  subplots  shows  the  cell  voltages  along  with  the  nominal  and  apparent  current 
densities.  The  24  individual  cell  voltages  are  displayed  in  thin  faint  colored 
lines  with  the  average  cell  voltage  based  on  the  measurements  in  a  thick  solid 
blue  line  and  the  estimated  cell  voltage  in  a  thick  dotted  red  line.  The  second 
row  of  subplots  show  the  oxygen  and  hydrogen  partial  pressures  in  each  GDL 
Section  (1-3)  as  well  as  in  the  channel  (4).  A  load  change  from  0.25  Acm_2to 
0.3  A  cm-2 occurs  at  approximately  183.4  min  followed  by  an  anode  purging 
event  at  approximately  183.7  min 

water  vapor  in  the  GDL  section  allowing  the  water  vapor  partial 
pressure  to  decrease  before  all  of  the  liquid  water  was  removed 
from  the  GDL  sections. 

The  liquid  water  saturation  in  the  GDL  section  closest  to  the 
channel,  san(3),  decreased  most  significantly  during  a  purge.  Liq¬ 
uid  water  flowed  from  the  GDL  towards  the  channel  until  the 
immobile  saturation  limit  was  reached,  san(3)  <  S{m,  at  which 
point  only  water  vapor  entered  the  channel  from  the  GDL.  Liq¬ 
uid  water  does  not  flow  from  the  GDL  to  the  anode  channel, 
following  the  purge,  until  the  liquid  water  saturation  in  the  GDL 
exceeded  the  immobile  saturation  limit.  If  the  purge  event  were 
to  have  occurred  over  a  longer  time  interval,  more  water  vapor 
in  the  anode  GDL  would  have  been  removed,  causing  a  more 
significant  impact  on  the  cathode  liquid  water  saturation  due  to 
the  water  vapor  transport  through  the  membrane. 
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Fig.  8.  Water  dynamics  during  a  load  change  and  purging  event.  The  first  row  of 
subplots  shows  the  water  vapor  partial  pressures  in  the  GDL  and  channels.  The 
channel  is  indicated  by  a  solid  line  and  the  three  GDL  sections  are  represented 
by  dashed  lines.  The  second  row  of  subplots  displays  the  liquid  water  saturation 
in  the  GDL.  Finally,  the  third  row  of  subplots  indicates  first  the  state  of  the  purge 
solenoid  valve  (0  indicates  the  valve  is  closed  and  1  means  the  valve  is  open), 
followed  by  the  mass  of  liquid  water  accumulating  in  the  anode  channel. 


8.  Model  validation  and  discussion 
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For  the  purposes  of  model  validation,  the  calibrated  model 
was  simulated  with  experimental  inputs  that  were  not  consid¬ 
ered  in  the  calibration  process.  The  resulting  model  predictions 
are  shown  in  Fig.  9  and  compared  with  the  actual  cell  volt¬ 
age  measurements  at  five  different  load  levels.  The  data  shown 
demonstrates  the  model  predicting  capability  over  a  range  of 
current  densities  and  air  stoichiometries.  At  approximately  162 
min,  the  air  stoichiometry  was  increased  from  200%  to  300%, 
causing  a  more  significant  increase  in  the  voltage  estimation 
(through  the  partial  pressure  of  oxygen  at  the  membrane  bound¬ 
ary)  than  was  experimentally  observed.  Despite  the  increased 
error  associated  with  the  oxygen  partial  pressure,  the  model 
correctly  estimated  the  degree  of  anode  flooding  at  various  cur¬ 
rent  densities,  correctly  predicting  no  significant  flooding  at 
low  loads.  As  the  load  level  was  reduced,  the  degree  of  flood¬ 
ing  decreased,  which  is  seen  from  inspection  of  the  difference 
between  the  apparent  and  nominal  current  densities  at  each  load 
level.  As  a  result,  the  deviation  in  voltage  decreased  between 
purges,  which  was  experimentally  confirmed. 

For  the  entire  validation  data  set,  the  average  estimation 
error  was  8.7  mV,  the  maximum  estimation  error  was  105  mV 


Fig.  9.  Model  validation  results.  The  first  subplot  shows  the  24  individual  cell 
voltages  in  thin  faint  colored  lines  with  the  average  cell  voltage  based  on  the 
measurements  in  a  thick  solid  blue  line  and  the  estimated  cell  voltage  in  a  thick 
dotted  red  line.  The  second  subplot  shows  the  nominal  and  apparent  current 
densities.  The  third  subplot  contains  the  anode  and  cathode  inlet  total  pressures. 
The  fourth  subplot  is  the  temperature  of  the  water  coolant  leaving  the  stack. 

and  the  standard  deviation  in  the  estimation  error  was  1 1.5  mV. 
Although  these  validation  error  statistics  are  approximately  two 
times  greater  than  the  error  statistics  associated  with  the  calibra¬ 
tion  data,  at  all  times  throughout  the  experiment  the  estimated 
average  cell  voltage  was  bounded  between  the  measured  mini¬ 
mum  and  maximum  cell  voltages  and  the  measured  cell  to  cell 
variation  was  larger  than  the  average  estimation  error. 

Although  the  model  of  the  reactant  and  water  dynamics 
results  in  an  accurate  estimation  of  the  voltage  degradation 
between  purges,  we  have  made  the  assumption  that  this  degra¬ 
dation  was  solely  due  to  the  accumulation  of  liquid  water  in 
the  gas  channels.  However,  it  is  conceivable  that  some  of  this 
degradation  could  be  due  to  the  accumulation  of  nitrogen  on  the 
anode  as  a  result  of  operation  with  air,  rather  than  pure  oxygen, 
or  catalyst  flooding.  Our  model  has  tunable  parameters  that  can 
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compensate  for  these  model  assumptions  and  simplifications, 
but  it  is  very  important  to  check  the  tuned  parameter  values 
against  other  published  values.  As  Table  2  shows,  the  tuned  aw 
and  tw i,  are  reasonable  and  within  the  range  of  published  results 
[6,27]. 

9.  Conclusions 

A  two-phase  isothermal  one-dimensional  model  of  reactant 
and  water  dynamics  has  been  developed  and  validated  using  a 
multi-cell  stack.  The  lumped  parameter  model  depends  on  six 
tunable  parameters  associated  with  the  estimation  of  voltage, 
the  membrane  water  vapor  transport  and  the  accumulation  of 
liquid  water  in  the  gas  channels.  During  step  changes  in  load, 
a  good  voltage  prediction  is  achieved  by  reproducing  both  the 
steady- state  and  dynamic  voltage  response  due  to  the  instanta¬ 
neous  increase  in  current  as  well  as  the  excursions  in  oxygen 
partial  pressure  resulting  from  the  manifold  filling  dynamics,  as 
demonstrated.  Finally  the  model  predicted  the  dynamic  effect 
of  temperature  on  voltage  as  shown  during  the  temperature 
transient  from  50°  to  60° C.  Although  simple,  this  model  cap¬ 
tures  the  voltage  dynamics  observed  in  a  fuel  cell  stack  at  low 
and  moderate  current  densities  under  the  range  of  conditions 
tested.  However,  caution  should  be  used  in  extending  this  model 
to  conditions  not  examined,  as  the  intent  of  this  model  is  for 
control. 

The  notion  of  apparent  current  density  is  a  means  for  describ¬ 
ing  the  impact  of  water  accumulation  on  the  dynamic  voltage 
behavior  of  a  PEMFC.  Future  work  is  focused  on  extending 
and  validating  this  simple  GDL  model  at  higher  current  density, 
and  establishing  the  connection  between  the  liquid  water  mass 
accumulation  and  voltage  using  neutron  imaging  techniques. 
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Table  A.l  ( Continued ) 


Symbol 

Definition 

AH  =  —228,  740Jmol-1 

Enthalpy  difference  from  STP 
(water  in  vapor  phase) 

('loss  =  1  mAcm-2  [51] 

Loss  current  density 

&ca,in  =  11.3E— ' 7  ms 

Cathode  orifice  constant 

&ca,out  =  11.3E— 7ms 

Cathode  orifice  constant 

&an,in  =  9.34E— 7  ms 

Anode  orifice  constant 

£an,out  =  9.34E— 7  ms 

Anode  orifice  constant 

kan,rm  =  11.31E-7mS 

Return  manifold  orifice  constant 

K  =  2.55E-13m2  [41] 

Absolute  permeability 

Ki  =  1.17 

Tuned  voltage  parameter 

K2  =  4.44  pA 

Tuned  voltage  parameter 

K3  =  1.78 

Tuned  voltage  parameter 

Ka  =  3.27 

Tuned  voltage  parameter 

L  =  3 

Number  of  GDL  sections 

Mh2  =  0.002  kg  mol-1 

Hydrogen  molecular  weight 

Mq2  =  0.032  kg  mol-1 

Oxygen  molecular  weight 

M^2  =  0.028  kg  mol-1 

Nitrogen  molecular  weight 

Mu 2o  =  0.018  kg  mol-1 

Water  molecular  weight 

n ceils  =  24  cells 

Number  of  cells  in  stack 

po  =  1  atm 

Standard  state  pressure 

R  =  8.314  Jmol-1  K-1 

Universal  gas  constant 

Sim  =  0.1  [41] 

Immobile  saturation 

AS  =  —44.43  Jmol-1  K-1 

Entropy  difference  from  STP 
(water  in  vapor  phase) 

tgd\  =  0.5  mm 

Total  GDL  thickness 

W>  =  0.038  mm 

PEMFC  membrane  thickness 
(includes  catalyst  layer) 

tw\  =  0.131  mm 

Tunable  water  layer  thickness 

T0  =  298.15  K 

Standard  state  temperature 

Vp  =  2.5  cm3 

GDL  section  pore  volume 

Vca  =  380  cm3 

Cathode  channel  volume 

Van  =  430  cm3 

Anode  channel  volume 

Van,rm  —  345  Cm3 

Anode  return  manifold  volume 

Qfyy  ==  1  1.5 

Tuned  water  vapor  diffusion 
parameter 

8y  =  0.167  mm 

GDL  discretization  width 

y  =  900  s— 1  [41] 

Volumetric  condensation  coeff. 

6  =  0.5  [41] 

Material  porosity 

<9C  =60°  [41] 

Contact  angle 

li  =  0.405  gm-U-1  [41] 

Liquid  water  viscosity 

p  =  997  kg  m-3 

Liquid  water  density 

cr  =  0.0644  Nm_1  [41] 

Surface  tension 

Anb,dry  —  1900  kg  m  3 

mb  dry  density 

37m b, dry  =  1.0kg  mol-1 

mb  dry  equivalent  weight 
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The  model  parameters  are  listed  in  Table  A.l  .  It  is  important 
to  note  that  for  convenience  the  parameter  values  listed  here  may 
have  units  other  than  those  expressed  in  the  modeling  equations. 


Table  A.l 

Parameter  symbols,  definitions  and  values 


Symbol  Definition 


Afc  =  0.030  m2 
bn  =  0.005 139[30] 
bn  =  0.00326[30] 

Dh2  =  114  mm2  s"1  [41] 
£>o2  =  30.3  mm2  s-1  [41] 
Ec  =  eekJmol-1  [51] 

F  =  96485  C  (mol 


Fuel  cell  nominal  active  area 
Ohmic  resistance  parameter 
Ohmic  resistance  parameter 
Hydrogen  diffusion  coefficient 
Oxygen  diffusion  coefficient 
Activation  energy 
Faraday’s  constant 
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